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A numerically stable method to solve the discretized Boltzmann-Enskog equation describing the 
behavior of non ideal fluids under inhomogeneous conditions is presented. The algorithm employed 
uses a Lagrangian finite-difference scheme for the treatment of the convective term and a forcing 
term to account for the molecular repulsion together with a Bhatnagar-Gross-Krook relaxation 
term. In order to eliminate the spurious currents induced by the numerical discretization procedure, 
we use a trapezoidal rule for the time integration together with a version of the two-distribution 
method of He et al. (J. Comp. Phys 152, 642 (1999)). Numerical tests show that, in the case of 
one component fluid in the presence of a spherical potential well, the proposed method reduces the 
numerical error by several orders of magnitude. We conduct another test by considering the flow 
of a two component fluid in a channel with a bottleneck and provide information about the density 
and velocity field in this structured geometry. 



I. INTRODUCTION 

Liquids often appear as homogeneous on a macroscopic scale, but not when observed on a microscopic scale where 
they may display density oscillations extending over a few molecular diameters. Equilibrium statistical mechanics 
theories such as density-functional theory (DFT) or integral equations can deal routinely with the presence of such 
inhomogeneities in density, concentration or other kinds of order parameters, and predict the ensemble average 
microscopic profiles and the associated surface and line tension, while a similar situation does not occur in non- 
equilibrium systems [THS]. In this case, the presence of inhomogeneities often causes difficulties in the numerical 
solution of the associated evolution equations. 

It is well known that the conventional hydrodynamic description, based on the Navier-Stokes equation, faces dif- 
ficulties when fluids are confined within a small volume or when the boundaries of the container have complicated 
shapes with typical lengths of the order of a few molecular diameters. Such a picture, while valid on a macroscopic 
scale, fails to describe very small systems [4H6]. On the other hand, the kinetic approach based on the distribution 
functions formalism and on the Boltzmann equation and its refinements represents a convenient description of both 
homogeneous and inhomogeneous systems. Among the existing numerical approaches employed to solve the Boltz- 
mann equation, the Lattice Boltzmann (LB) method plays a prominent role [7-9J. It is a discretized version of the 
continuous Boltzmann equation and gives good results in the homogeneous phases jlDHH]. However, the numerical 
solution of inhomogeneous systems within the LB scheme is challenging: as reported by several authors p!3HT8] , a 
straightforward application of the LB equation leads to the observation of an unphysical effect, the so-called spurious 
currents, resulting from the discretization procedure. To cure this pathology, molecular interactions must be handled 
with care. 

In the literature, internal forces are accounted for in two different ways: a) either by imposing the condition that 
the equilibrium distribution gives the desired form of the pressure tensor or b) by introducing an appropriate forcing 
term [19', (20]. The forcing term can be chosen in two different manners, either proportional to the gradient of the 
pressure excess over the ideal gas value or proportional to the product of the density times the gradient of the excess 
chemical potential, that is, by using the Gibbs-Duhem condition in differential form. Actually, the second choice is 
consistent with microscopic theories, such as DFT [Ij, where the equilibrium condition is given by requiring that the 
gradient of the local chemical potential is locally balanced by the external forces. 

In the present paper, we discuss a LBE algorithm based on the Boltzmann-Enskog transport equation |2TH25] . The 
approach is particularly convenient when the packing effects are relevant, that is from moderate to high fluid densities. 
A straightforward application of the LBE algorithm leads to numerical instabilities so we introduce a numerical scheme 
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that employs a trapezoidal time discretization plus an extension of a procedure, originally proposed by He et al [^B], 
that uses two distribution functions, instead of one, to reduce the spurious currents phenomenon. In this scheme, 
one distribution function tracks the local density profile while the other tracks the local momentum density. The 
standard phase space distribution function /(r, v,t) is evolved concurrently with an auxiliary distribution function, 
named ^(r, v,t), whose zeroth velocity- moment is the hydrodynamic pressure and its first moment is identical to the 
corresponding moment of /(r, v, t). According to previous authors, the reason for the increased stability of the double 
distribution method stems from the fact that the forcing term in the g-equation is multiplied by the difference between 
the local and the global Maxwellian thus reducing its importance with respect to the original f-equation, where the 
forcing term is multiplied by the local Maxwellian. The method was later extended and generalized by T. Lee and 
coworkers [27H29] . 

The main difference between our approach and the previous ones, besides the bottom-up microscopic modeling of 
the fluid proposed in earlier work [32 - 34J, consists in the choice of the function employed to define the g-distribution 
function. As we shall see, with the present choice it is straightforward to generalize the method to multicomponent 
fluids, while in the original formulation such a generalization is not straightforward. In this way, our method leads 
naturally to a form of the forcing term similar to that in the Gibbs-Duhem route. This strategy can also be generalized 
to multicomponent fluids, whereas the pressure route cannot. 

The paper is organized as follows: in Sec. |TT]we present the evolution equation for the one particle distribution 
function / and for the auxiliary distribution g both for the simple fluid and for the fluid mixture. In Sec. Ill we 
discuss the discretization procedure In Sec. |IV| we present numerical tests of the proposed method. Finally in Sec. |V| 
we present our conclusions and perspectives. 



II. EQUATIONS FOR THE DOUBLE DISTRIBUTION FUNCTIONS 

We start the discussion with the set of Boltzmann-Enskog equations characterizing a mixture of M species, labelled 
with an upper index a = 1, M. The evolution equation for a particular distribution function /"(r, v, t) can be written 
as: 

gr(r,v,t) = ■ ^r{r,v,t) + ^ J"/^(r,v,t) 

(1) 

where the material derivative is given by: 

and F"(r) is an external velocity independent force field acting on component a and J"^, represents the effect on the 
single particle distribution function of the interactions among the fluid particles of type a and p. 

Using a separation of the interaction term into a kinetic rapidly varying part and an hydrodynamic part originally 
introduced by Santos et al. 122 and extended to mixtures later ^51 [36] we rewrite ([T]) as: 

^ = -ujir-f^^) + Sf{v,v,t) (3) 

The first term in the r.h.s. of eq. ([3| is a Bhatnagar-Gross-Krook (BGK) relaxation term [37], uj an inverse 
relaxation time, and Sj is a source term due to external forcing and molecular interactions. According to ^SSJ it can 
be written as 

Sf (r, V, t) = - 5^ A/« + ;3(v - u) ■ (r, t)r„(r, v, t) (4) 

with 13 = l/Zc^T, T is the temperature and the Boltzmann constant. In addition, is a Maxwellian velocity 
distribution whose mean velocity is the local fluid velocity u(r,t): 

where mv^ = ksT for particles of common mass m, and 

f:^{v, V, t) = n"(r, V, + /?[K(r, t) - u(r, t)) ■ (v - u(r, t))] }r„(r, v, t) (6) 
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with u*^ the average velocity of the component a. The term C*^ is a cohisional kernel describing the change of due 
to the interactions. 

We first rewrite ([3| in a form that is equivalent up to terms of third order in the Hermite expansion 

S'} = I3 (^C" + n"?^^ • (v - u)r„. (7) 
From the phase space distribution function /"(r,v,t) we can compute the particle partial density, 

n"(r,t) = j dvr(r,v,t) (8) 
and the momentum current carried by particles of type a, 

n-(r,t)u-(r,t) = j dvr(r,v,t)v (9) 
which from eq. (|3| satisfies the continuity equation 



The average fluid velocity is obtained from 



^+V-Ku") = 0. (10) 
ot 



with the global density given by 

The numerical solution of eq. ([3| is plagued by numerical instabilities as reported in ref. [26], because the term SJ 
featured in the r.h.s. is quite large in the interfacial regions, since the main contribution to C^, which is proportional 
to the gradient of the local chemical potential, varies rapidly. Alternatively, following the seminal idea put forward 
by He and coworkers in ref. [26j and pursued by Lee and coworkers [29j to stabilize the numerical solution the one 
component version of eq. ([3|, it is possible to employ an auxiliary distribution function, named ^"(r, v,t) such that 
the role of the forcing term featured in its evolution equation is effectively reduced. Such an heuristic recipe stabilizes 
the numerical solution by "decoupling" the density and the momentum equations. In the present treatment, we will 
handle the stabilizing terms in an effective way, without relying on any heuristics. Let us introduce the auxiliary 
distribution 

^-(r, v,t) = r(r, v,t) + (n-(r,t) - n^{v,t))V^ (11) 

where n"(r, t) is a function of the partial densities , to be determined in the following, and Fq indicates the velocity 
distribution at global equilibrium, that is, the Maxwellian corresponding to u = 0. One assumes that the function 
depends from its argument through {n"(r,t)}. From the definition (11) one can see that differs from w.r.t. 
the zeroth moment 



^ dv5"(r,v,t) = n"(r,t), (12) 

but shares the same first moment 

''dvc/"(r,v,t)v = n"(r,t)u"(r,t). (13) 



By using eqs.(ll) and ([3|, the evolution equation for ^"(r, v,t) reads 

Dg"^ D/" D 



Dt Dt Dt 



(n--n")Fo. (14) 



— (n"(r,t) - n«(r,t)) = (v - u) • V(n« - n«) - ^„''(— - 5^p)V ■ - ^(^ - 5^p){u^ - u) • Vn" (15) 

/3 /3 
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one obtains the evolution equation (14) for g^{T,v^t) as 



with 



5,"(r,t)=/3(c" 



Dt 



.n"— .(v-u)(r, 

m 



(16) 



m J 



•(v 



u)ro-^n''( 

13 



and 



g^^{r, V, = n"(r, t) (l + /?m(u"(r, i) - u(r, t j) ■ (v - u(r, t)))r„ + (n«(r, t) - n"(r, i))ro. 



(17) 



(18) 



It can be checked that the evolution equation for g'^ or the one for lead to the same balance equation for u". In 
practice, in the numerical work we shall use the equation to determine the density and track the formation of 
interfaces, and the g^ equation to determine the velocity field u". 

The main motivation behind the transformation from to g'^ is that the effect of the forcing term, 5*^, featuring 
in ([l6| can be rendered smaller than the corresponding effect due to the forcing term, SJ, in the original equation 
([3| for by an appropriate choice of the function n"(r,t). In fact, the first term in Sg is of order (v — u)^ because 
it contains the product of (v — u)(r^^ — Fq) , whereas the second term can be rendered small using the arbitrariness 
of the function 11" is. As far as the last term is concerned we shall verify that the last term in Sg is actually small in 
our numerical simulation. One expects that a weaker forcing term helps the stability of the numerical solution. In the 
one component case He and coworkers suggested to replace /3~^n" by the thermodynamic pressure pt. In order to see 
that we use the explicit representation of the function C", which represents the effect of the molecular interactions 
in the model studied. 

We first separate the effective field into three separate contributions, the separation being quite generic and not 
determined by the particular model used: 



The first term can be written as: 



C".™/ = _„"(r,i)V<„,(r,i). 



(19) 



(20) 



where /if^^ is the non ideal part of the chemical potential of the a component. For density profiles smooth enough we 
can write 



/3 



and for the viscous part 



Cr''^{r,t) « -n"(r,t) ^(r/^'^V^wf (r,t) + {vf + l.v"^)V,{V ■ u^)) 



(21) 



(22) 



The coefficients 7 and r] depend on the specific model. In appendix A, we report their explicit representation for a 
system of hard-spheres with attractive interactions. 

In the case of a one component fluid it is straightforward to derive the equation for the g distribution which closely 
resembles the equation derived by He and coworkers. After dropping the unnecessary index a one has: 



Dg 

Dt 



= -ojig - geq) + Sg 



(23) 



Sg = /3{C + n-) ■ {v - u){T, 

m 



To)- 



Vn - Vn + ;5C ■ 



/3F' 



(v - u)ro - n( 



da 

dn 



i)(v-u)ro. 



(24) 
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Using (20) and neglecting the non equilibrium contributions to C we have: 

- Vn + /3C = -/3nV/i (25) 

where /i is the total chemical potential. Finally with the help of the Gibbs-Duhem relation we introduce the thermo- 
dynamic pressure : 

Vpt = nVii. (26) 



Hence, requiring the vanishing of the square brackets in (24) is equivalent to the condition: 



/3Vn = Vpt-n — . (27) 
m 

In other words choosing j3Ii to be the thermodynamic potential augmented by the contribution due to the external 
field makes the second term of ( [24| ) to vanish. From the physical point of view such a condition is a consequence of 
the hydrostatic equilibrium condition [Ij. 

Unfortunately, in the multi-component fluid the identification of with the pressure is not possible. The reason 
is that in this case the distribution functions one needs a 11" function for each component, whereas one can find only 
one pressure, through the Gibbs-Duhem relation 

Vp, = ^n"V/i". (28) 

a. 

Moreover, by using the pressure route it is very difficult to obtain a satisfactory numerical solution in the general 
case, as in presence of confining walls, spontaneous layering mechanisms, or free interfaces. Alternatively, we choose 
the unknown function 11" as the "potential function" associated with the vector field Vn" — /3C", in such a way as 



to cancel this term from eq. (17). More precisely, the function 11" is chosen to be: 

n"(r,t) = n"(r,t)-/3 j dr' (c"'^'^(r^ t) + n"(r^ t) (29) 
Therefore, being C"'^-^ a functional of density, U" is chosen to be a non-local function of density, in stark contrast 



with previous proposed approaches that are based on a local compensating pressure term [26| [291431] . Eq. (29) also 
provides the operational route to our approach. In fact, the integral is evaluated numerically using trapezoidal spatial 
integration, which provides a satisfactory numerical solution in terms of accuracy. It should be noticed that, being 
an integral over a vector field, the integration depends on the origin and the specific path of the integral. However, 
this aspect is not dangerous for systems where a symmetry point can be found. In addition, the integration constant 



never appears in the evolution equation and thus does not need to be determined. Although eq. (16) looks more 
complicated than the original one, it behaves better in numerical terms and gives rise to smaller interfacial currents, 
as shown in the sequel. 

III. NUMERICAL SOLUTION 

We illustrate the numerical solution of the proposed method by considering explicitly the one component case, 
while the multicomponent case can be easily deduced. Let us consider again the integration of the generic evolution 
equation 

^=l^(/,M)(r,v,t) (30) 

where the unspecified kernel contains both the collisional term, the BGK term and the external force ¥ • d^f. As 
customary in the derivation of the Lattice Boltzmann method, the distribution function is first projected on an finite 



Hermite basis set to handle the dependence on velocity [39| |4Q| . By taking eq. (30) as our reference equation, the 
r.h.s. depends on / but also on its moments M = {M^}, with 

Mp{r,t)=< f\np> (31) 

where Hp is the p-th Hermite polynomial, and 

<A\'Hp>= J dvA{r,v,t)Hp{v) (32) 
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expresses the Hermite scalar product. 

In order to discretize eq. (30), we start by considering the following truncated Hermite expansion 

K 



/(r,v,t) =ro(«)X^^^Mj,(r,t)Hp(v) (33) 



P- 



where K is the order of truncation of the Hermite expansion and Mp =< f\Hp >=< f\l-Lp >. In fact, from the 
definitions, it follows that the original and the truncated forms of the singlet distribution share the same moments up 
to p < K. By the same token, we consider the expansion of the collisional kernel. 



f^(r, V, t) = Toiv) f2 ^0,(r, t)n^^^ (v) (34) 



with Op =< Q\H^P'' >. As for the distribution function, Q has moments O = {Op} shared by the full and truncated 
representations of the kernel ft. 

The LBM is based on replacing the Hermite scalar products by Gauss-Hermite quadratures to evaluate its moments, 

G 

Mp =< f\Hp >=Y. fp^^""^ (^p) (35) 

where the vectors Cp are a set of quadratures nodes, Wp are the associated weights, and G is the order of the 
quadratures. 

The operational version of the LBM scheme is provided by the following quantities /p(r, t) = Wpf{T,Cp^t)/ro{cp) 
and Qp{r^t) = WpQ{r^Cp^t)/To{cp). From these transformations, the evolution equation of the new representation 
reads 

^/j,(v,t) + Cp • V/p(v,t) = np{fp,M){r,t) (36) 

where we have rewritten the streaming term v • V/ in its Hermite form. The exact time evolution of the populations 
over a timestep h then reads 

rt+h 

fp{r^Cph,t^h) = fp{r,t)^ J ds%{fp,M){v,s) (37) 
On the other hand, a second-order accurate 0{h^) numerical integration can be obtained via the trapezoidal rule [41j, 
dsQpifp, M)(r, s) = - [Qpifp, M)(r + Cph, t ^ h) ^ %{fp, M)(r, t)] + 0{h^) = - (1^^+^ + ^1) + 0{h^) (38) 



where = Vt{fp,M){v,t). 

Eq. (|38| is apparently implicit. However, the scheme can be rendered explicit by using the following exact mapping 
to transform the original populations into the new set 

The moments in the / representation, collectively called Mp = {< f\l-Lp >}, are related to those in the / represen- 
tation by 

Mp = Mp - (40) 



In many circumstances, both relations eqs. (39)-(40) are invertible, that is, we can obtain explicitly the populations fp 
as a function of fp and the moments M as a function of M. This is the case, for example, for BGK or Fokker-Planck 
kernels [42 J, and in presence of external forces. 

Finally, the temporal evolution for the populations fp if given by the following updating scheme 

fp{v + Cph, t^h)= fp{r, t) + %{fp, M)(r, t)h (41) 
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that provides the way to integrate the equation via the trapezoidal route. It is practical to work in the / representation, 



and substitute the quantities fp and M in the collisional kernel featuring in the r.h.s. of eq. (41). 

For the collisional kernel C appearing in eq. ([3|, however, the relation is non- invert ible, since C is a functional of 
the hydrodynamic moments. Yet, by decomposing 0{M) = 0^^^{M) + C[M], where O^^^(M) is the residual part 
of the collisional moments, being a function (rather than a functional) of the hydrodynamic moments, a workable 
algorithm is obtained via the scheme 

^,=„-m._aM, ,42) 



so that the original moments M are expressed as functionals of M, and substituted in eq. (41 ). It is straightforward to 
show that, for C = 0, the BGK and external forcing components give rise to the second- accurate integration method 
introduced by Guo et al. [42j . 



IV. RESULTS 



In the following, we will analyze an ideal fluid with both the Euler integration, referred to as EU, and trapezoidal 
integration, referred to as TR. Subsequently, we will consider the hard-sphere system and compare the simulations 
obtained via the single distribution method (without the auxiliary distribution), named SD, and the double distribution 
method, named DD. By distinguishing the case of Euler integration from the trapezoidal integration, we have four 
combinations, for example, the double distribution method with the trapezoidal rule will be named DD-TR, and 
analogously for the other combinations, such as SD-EU, SD-TR and DD-EU. 



A. Ideal fluid in a potential well 



To illustrate the numerical capabilities of the LBM, let us first consider an ideal fluid (by setting the collisional 
kernel C = 0) in the presence of an external central potential, expressed as 



ext 



I -6(1 



+ cos(^)) ifr<^ 



(43) 



where the potential depth is taken to be e = 0.3 x ksT and the well size ^ is varied in order to compare the standard 
Euler versus the trapezoidal integration rules. The external force is expressed as F^^^ = — acting on particles 

of unit mass. At global equilibrium, the density should be distributed as n^^ ex.p{—U^^^ /v^) , with no = y /^(irn(r), 
and the current should be zero everywhere. 

By applying the trapezoidal rule, the populations fp are updated in time and, at every time step the density and 
current are computed as n = ^p fp and J = ^pCpfp. In the / representation, the hydrodynamic moments that 
contain the second-order accuracy in space and time are computed by reversing equation (42), so that n = fp = n 

and J = EpCp/p=J + F^-*|. 

Finally, the expression of the external forces up to second Hermite order, reads 



pext . ^(1) _^ 2F^^*u : -HI^) 



(44) 



where ni^^ = H^^\cp) = ^ and H),"^ = n^^\c^) 
and I is the unit tensor. 



By using eq. ^ it follows that /j, = (1 + x [/p + ^f^i 

following post-collisional term 



, being a vector and a tensor of rank two, respectively, 
Tp] and the populations are updated to the 

h 



1 



to be contrasted with the standard Euler integration, reading 

/; = fp + cvh [f;i{n, u) - fp] + hj^; 



-pext 

ujhl2 P 



(45) 



(46) 



Both the trapezoidal and Euler evolutions are then completed by the streaming stage, reading /p(r + hcp^ t -\- h) 
fp{v^ t) and /p(r + hcp^ t -\- h) = fpiv^ t), respectively. 



t/Ax t/Ax 

ABC 

Figure 1: Numerical error in the density (panel A), velocity (panel B) and current (panel C) for the ideal gas system in presence 
of the central potential well. See text for details. The dashed and dot-dashed lines represent the power law dependence of the 
numerical error, as reported in the legends. 



We simulate a three-dimensional system and in Fig. [T]we report the error on density as Err(n) = maxr(|n— n^'^|/no), 
the error on fluid velocity arising from parasitic effects, as Err{\x) = maxr(|u|/vT), and the error on current as 
Err (J) = maxrd J|/no^'T)- The data show that the numerical errors in the density, velocity and current decrease 
systematically with the mesh resolution for both the EU and TR methodologies. The error is reduced by about 
two orders of magnitude for the TR method as compared to the EU scheme. In particular, the error in density 
decreases as Ax^ for both methods, while the error in current drops as Ax^ and Ax^ for the EU and TR methods, 
respectively. These preliminary results provide a reference for the subsequent simulations of the hard-sphere system 
and an important indication on the quality of the trapezoidal evolution method. 



B. Hard sphere fluid mixture in a potential well 

We now consider a non ideal fluid mixture of hard spheres, and numerically solve the statics of the problem in the 



presence of the same central external potential of eq. (43), and integrate the dynamics with and without the auxiliary 
distribution method. 

The trapezoidal integration for the two distributions generalizes eq. ( [41] ) to 

f^(r + Cph,t + h) = f^ir,t) + hni^if^,{M})ir,t) (47) 

g;{r + Cph,t + h) = g;{r,t) + hni^{g;,{N}){r,t) (48) 

where 

^hifp, {M}) = - /;) + sip (49) 

^gA9p, m) = ^idp'^' - 9p) + S% (50) 

Here, {M} and {TV} refer to the set of moments of the populations and respectively. 
We compute the relevant moments, that is, densities, HS chemical potentials and currents as 

n" = Y.fp (51) 

P 

n" = Y.9p (52) 



J" = nfi" = ^cj,5^ (53) 

P 
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Figure 2: Numerical error in the fluid velocity in the presence of the central potential well for the one component fluid (left 
panel) and the binary mixture with diameters of cta = 4 and as = 8 (right panel). 



then 



p 



P 



The explicit form of the r.h.s. of eqs. (49)- (50) reads 



a,eq 



a,eq 



9jP 



" + n"u" • n^^^ + n"(2u"u - uu) : n^^^ 
Wp + n"u" • n^^^ + n"(2u"u - uu) : n^^^ 



(54) 
(55) 
(56) 

(57) 
(58) 
(59) 

(60) 



In Fig. |2j the numerical error in the computed velocity profiles is reported for the one component and the two 
component fluids. The error decreases with increasing resolution and is smaller by respectively a factor 10 and 50 
for the case of SD-TR and DD-TR simulations as compared to the SD-EU method. The data are similar for the one 
and two-component systems, follow the same behavior observed for the ideal gas, and the error in velocity decreases 
steadily with increasing resolution. The spurious velocities are about 50 times smaller for the DD-TR case as compared 
to the SD-EU integration. 

A major advantage of the trapezoidal integration alone is the possibility to work at high packing fractions, up 
to about 0.35, whereas with standard Euler integration, the maximal packing fraction before numerical instabilities 
develop is 0.27. 
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Figure 3: Sketch of the channel flow system in presence of an obstacle. The channel is filled with either a one- component fluid 
or a binary mixture. For the simulations we have set the geometry to L = 120, / = 30, H — 35 and /i = 28 in lattice units. The 
diameters of the large and small particles are denoted a a and as, respectively. 



C. Channel flow with a bottleneck 

We now consider the flow of a one-component HS fluid and a binary mixture in a channel flow, in the presence of a 
bottleneck, as depicted in Fig. [3j Flows in the presence of a sharp obstacle represent a critical test to the numerical 
methodology due to the harsh collisions that the particles experience with the corners of the obstacle. In particular, we 
choose a rather strong forcing term, being equal to 10~^ in lattice units, in order to obtain large impinging velocities 
against the obstacle. We further impose no-slip boundary conditions on the fluid populations at the solid wall for 
both the / and the g distributions. For this we employ the mid-point bounce-back rule on the populations [7J. We 
initially simulate a one-component system composed of hard spheres of diameter <j = 8 and make complementary 
simulations with a two-component mixture with hard spheres diameters of = 4 and = 8. 

As Fig. |4] demonstrates, the naive SD-EU method provides strong spurious velocities arising from the presence of the 
wall. In fact, away from the obstacle, the streamlines are expected to be parallel to the wall, whereas we observe strong 
non-parallel streamlines near the wall that confirm the low quality of this type of simulations. Conversely, the DD-TR 
method provides well aligned streamlines near the wall and far away from the bottleneck. From these observations, 
we decided to consider further benchmarks by looking at the results obtained with the DD-TR methodology alone. 

Hard spheres in proximity to an irregular surface present an interesting phenomenology in itself. In fact, in proximity 
to the obstacle, the fluid particles go around the obstacle with non-trivial patterns. In particular, as the flow lines 
in Fig. [5] reveal, a first bounce back is found near the convex corner. Entropic forces have a strong influence on the 
spatial distribution of the particles and, as previous studies demonstrate ^43j, the concave corners effectively attract 
particles, while convex corners exert repulsive forces. Such dual behavior is recovered by our simulations, as revealed 



11 




120 



1-1 1.0 



0.8 



0.6 



0.4 



120 



r 

■o.o 



Figure 4: (Color online) One component system: comparison of streamlines for packing fractions of 0.13 as obtained with the 
SD-EU (upper panel) and with the DD-TR (lower panel) methods. 



by the density profiles in Fig. [5] where the accumulation of particles toward the edges of the obstacle is clearly visible. 
In addition, we observe that the density profiles have a very weak dependence on the forcing term, with a somehow 
stronger variation in proximity of the corners for the incoming particles, as compared to the static case. A further 
validation of the method is given by the computation of the divergence of velocity, as reported in Fig. [6j For the 
compressible system considered here, the quantity V • u should be zero everywhere, while spurious compressibility 
effects are clearly visible when employing the SD-EU method. The DD-TR method minimizes such error up to three 
order of magnitudes. 

For the system at hand, the simulations provide new interesting information about the fluid velocity in this geometry, 
as shown in the following. 

The simulations provide the fine details of flow pattern for the one-component system at varying packing fraction, as 
illustrated in Fig. [7[ For increasing packing fraction, the streamlines become more and more disordered in proximity 
to the convex corners of the bottleneck. A quite disordered pattern is observed already at a packing fraction of 
0.26, with flow separation appearing in correspondence with the impinged corner. The dynamical disorder appears 
to initiate at the far away edge of the obstacle with respect to the incoming flow direction. At a packing fraction of 
0.34, the disorder has propagated to the whole region of the bottleneck with rough recirculation patterns. It should 
be noticed that, for increasing packing fraction, the modulus of velocity is reduced overall, with strong peaks localized 
near the corners. 

We have next considered a binary mixture of hard spheres of diameters cja and (Jb-, flowing in the same channel 
with the bottleneck. An important aspect of the binary mixture is that entropic forces play different roles on the 
species with different diameters. For particles of smaller diameter, entropic forces are smaller, and these particles can 
distribute more uniformly between the concave and convex corners. Consequently, the flow pattern is expected to be 
more ordered. This behavior is shown in Fig. [Sj for a binary mixture with = 8 and = 4. The streamlines of 
both the large and small particles have a smoother behavior as compared to the one-component case. The modulus 
of velocity of both species is more uniform as compared to the one-component system, with a smoother distribution 
around the obstacle. In Fig. |9) the binary mixture with particles of size cr^ = 8 and (Jb = 2 presents even smoother 
flow lines and smoother distribution of the velocity moduli, as compared to simulations at smaller size ratio and at 
the corresponding packing fractions. Overall, we conclude that in the binary mixture, the component with particles 
of smaller size acts as a powerful lubricant that regularizes the flow pattern and distributes evenly the flow velocity 
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Figure 5: One component system: density profiles for different values of the x coordinate as obtained with and without forcing. 
The profiles for x = 40 and x = 80 correspond to positions right before and right after the corners of the obstacle. For all x 
values, the profiles with and without forcing are basically indistinguishable. Profiles have been shifted upward for the sake of 
clarity. 




10 20 



z 



Figure 6: One component system: divergence of velocity D = V - u/{vt/^x) computed at mid channel {x = 60) for the system 
in fiow condition. The two profiles correspond to the DD-TR (circle symbols) and SD-EU (square symbols) methods. The inset 
displays the ratio \d^d-tr ^ jjSd-eu 



over the whole system. 
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Figure 7: (Color online) One component system: streamlines for packing fractions of (A panel), 0.13 (B panel), 0.26 (C 
panel) and 0.34 (D panel). The color map represents the modulus of the flow velocity, normalized by its maximum. 



CONCLUSIONS 



In this paper, we have illustrated a numerical version of the Lattice Boltzmann method for the simulation of hard- 
sphere one-component and binary mixtures, that can deal with rapid spatial variations in the number density. As 
well-recognized in the Lattice Boltzmann community, strong inhomogeneities in the density induce strong parasitic 
currents that need to be handled with great care. 

In our method, we have extended the previous ideas of He et al. [26j and T. Lee [29j but with some important 
modifications. In particular, we compute the excess chemical potential arising from the hard sphere collisions on- 
the-fly, without resorting to an educated guess of its functional form. In addition, we have adapted the trapezoidal 
integration rule for the time evolution of the populations, written as an explicit time-stepping algorithm. 

The numerical results showed that, at all packing fractions considered in the benchmarks, the method provides 
robust results and stable numerical behavior. 

We conclude by mentioning that the presented method can be applied without major modifications to nanofluids in 
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Figure 8: (Color online) Binary mixture with aA — S (left column) and as = 4 (right column): streamlines for packing 
fractions of 0.0 (panel A), 0.13 (panel B), 0.26 (panel C) and 0.34 (panel D) and for 50% composition. 



presence of electrostatic interactions, as presented in ref. [33]. For these systems, internal electrostatic forces exerted 
between charged species arise from the solution of a Poisson problem treated at mean-field, Vlasov level. Also in this 
case, the trapezoidal and double distribution methodology can be applied straightforwardly since electrostatic forces 
are treated at the same level of external forces. 



VI. APPENDIX 



We report the formulae given elsewhere which have been used to compute the various terms of the effective field. 
The details have been reported in a previous publication [36j. In eq. (19) we can identify a force acting on the a 
particle at r due to the influence of all remaining particles in the system, the so called potential of mean force. For a 
hard-sphere mixture we have: 

C«,-/(r, t) = -kBTn^ir, t) ^ cj% j dkk^«^(r, r + a,^k, t)np(v + a^^k, t) + n-(r, t) ^ G-^(r, t) (61) 

with (JcK^ = (cTc^a + cr/3/3)/2, while the last term represents the molecular fields associated with the attractive forces: 

G"^(r,t) = - j dr'n^(v',t)g^^(v,v')VU''^(v-v') (62) 
with [/"^(r) a long range attractive potential. The drag term is: 



C"-^™9(r, t) ^ -n"(r, t) ^ ^al^^j'^ ^g^^({n-{r, t)})nP{r, i)(u"(r, t) - u^{r, t)) (63) 




and for the viscous part 



where g^fd is the pair correlation function evaluated at contact (r = dafd) As shown in Ref. [36] one can derive the 
following expressions in the limit of a uniform system for the viscosity: 



15 "^V TT 

and 



aB 47r 4 mkBT ^ 

= T^^a/3\/ -^^9ai3n^ (65) 



Vf = Iv"-^ (66) 
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